#include <lineages.h>
#include <vector>
#include <string>
#include <map>
#include <iostream>
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <algorithm>
#include <boost/random.hpp>
#include <boost/random/gamma_distribution.hpp>
#include <boost/random/mersenne_twister.hpp>

template <typename T>
inline double sum(vector<T> all)
{
    double res = 0.0;
    for(unsigned int i=0; i<all.size(); i++)
    {
        res += (double)all[i];
    }
    return res;
}

template <typename T>
inline double mean(vector<T> all)
{
    double res = 0.0;
    if(all.size() == 0) return res;
    for(unsigned int i=0; i<all.size(); i++)
    {
        res += (double)all[i];
    }
    return res/(double)all.size();
}


template <typename T>
inline double stdev(vector<T> all)
{
    double res = 0.0;
    if(all.size() == 0) return res;
    double avg = mean(all);
    for(unsigned int i=0; i<all.size(); i++)
    {
        res += pow((double)all[i] - avg,2);
    }
    return pow(res/((double)all.size() - 1),0.5);
}

inline double random(double start, double end)
{
    return start +(end-start)*rand()/(RAND_MAX + 1.0);
}

ostream& operator<<(ostream& os, const Lineages&)
{
    return os;
}

Lineages::Lineages(vector<int> terminalTypes,vector<string> terminalLineages)
{
    //int testInt;
    termTypes = terminalTypes;
    //testInt = termTypes.size();
    termLineages = terminalLineages;
    //testInt = terminalLineages.size();
    linCompl = -1.0;
    linRobust1 = -1.0;
    linRobust2 = -1.0;
    linRobust3 = -1.0;
    depthMax = -1;
    depthMean = 0.0;
    depthSD = 0.0;
    cellDeletionRate = -1.0;
    ruleDeletionRate = -1.0;
}

Lineages::~Lineages()
{
}

void divideTerminal(string termLineage,int childType1,int childType2)
{
}

double Lineages::complexity()
{
    if(termTypes.size() <= 1) return 0;
    if(linCompl>0.0) return linCompl;
    coalescent();
    return linCompl;
}

double Lineages::robustness(int disturbType, int robustType)
{
    if(disturbType == 1)
    {
        return robustness1(robustType);
    }
    else if(disturbType == 2)
    {
        return robustness2(robustType);
    }
    else if(disturbType == 3)
    {
        cerr << "not defined" << endl;
        exit(1);
        return robustness3(robustType);
    }
    else if(disturbType == 22) // use sampling to get robustness2
    {
        return robustness2Sample(robustType);
    }
    return 0.0;
}

double Lineages::robustnessGamma(int disturbType, int robustType, double shapeK, double meanU)
{
	if(disturbType == 1) // type 1 based on gamma distribution death rate
    {
		return robustness1Gamma(robustType,shapeK,meanU);
    }
    else if(disturbType == 22) // type 22 based on gamma distribution of fail rate
    {
		return robustness2SampleGamma(robustType,shapeK,meanU);
    }
    return 0.0;
}

double Lineages::robustness2SampleGamma(int robustType,double shapeK,double meanU)
{
	if(linRobust2 > 0.0) return linRobust2;
    if(internalTypeCnt <= 1) return 0.0;
    if(allLineages.size() < 1) coalescent();
    meanU = ruleDeletionRate;
    if(ruleDeletionRate < 0)
    {
        cerr << "Define ruleDeletionRate first!" << endl;
        exit(1);
    }

    int sampleCnt = allLineages.size()*10;

    boost::mt19937 rng;
    boost::gamma_distribution<> gd(shapeK);
	boost::variate_generator<boost::mt19937&,boost::gamma_distribution<> > var_gamma(rng, gd);

    vector<int> allRulesOri;
    for(int i=0; i<(internalTypeCnt-1); i++) allRulesOri.push_back(i);
    //
    // use (internalTypeCnt-1) so that you won't remove the rule governing the first division from the fertilized egg
    //

    vector<int> pickedRules;
    vector<int> pickedCells,tmpVec;
    double rndDbl;
    int failCnt = 0;

    map<int,vector<int> >::iterator it;
    map<int,int>::iterator it2;
    map<int,int> totalTermType2Count; // total number of each types at terminal
    vector<int> affectedDescendant;
    vector<int>::iterator it3;
    map<int,int> affectedTermType2Count; // number of each types at terminal of a sub-lineage
    vector<double> fitnessUponDeath;

    for(unsigned int i=0; i<termTypes.size(); i++)
    {
        if(termTypes[i] != sparableCellType)
        {
            it2 = totalTermType2Count.find(termTypes[i]);
            if(it2 == totalTermType2Count.end())
            {
                totalTermType2Count[termTypes[i]] = 0;
            }
            totalTermType2Count[termTypes[i]] ++;
        }
    }
    if(totalTermType2Count.size() < 1) // none of the termini is important
    {
        return 0.0;
    }

    for(int i=0; i<sampleCnt; i++)
    {
    	failCnt = 0;
        pickedRules.clear();
        affectedDescendant.clear();
        affectedTermType2Count.clear();
        for(unsigned int j=0; j<allRulesOri.size(); j++)
        {
        	rndDbl = random(0,1);
        	if(rndDbl < ruleDeletionRate / shapeK * var_gamma())
        	{
        		it3 = allRulesOri.begin() + (ptrdiff_t)j;
        		pickedRules.push_back(*it3);
        		failCnt ++;
        	}
        }

		//cerr << ((double)failCnt / (double)allRulesOri.size()) << endl;

        for(unsigned int j = 0; j<pickedRules.size(); j++)
        {
            pickedCells = internalTypes2allIndex[pickedRules[j]];
            for(unsigned int k = 0; k<pickedCells.size(); k++)
            {
                tmpVec = allIndex2DescendantIndex[pickedCells[k]];
                affectedDescendant.insert(affectedDescendant.end(),tmpVec.begin(),tmpVec.end());
            }
        }
        sort(affectedDescendant.begin(),affectedDescendant.end());
        it3 = unique(affectedDescendant.begin(),affectedDescendant.end());
        affectedDescendant.resize(it3 - affectedDescendant.begin());
        for(unsigned int j=0; j<affectedDescendant.size(); j++)
        {
            if(termTypes[affectedDescendant[j]] == sparableCellType) continue;
            it2 = affectedTermType2Count.find(termTypes[affectedDescendant[j]]);
            if(it2 == affectedTermType2Count.end())
            {
                affectedTermType2Count[termTypes[affectedDescendant[j]]] = 0;
            }
            affectedTermType2Count[termTypes[affectedDescendant[j]]] ++;
        }
        fitnessUponDeath.push_back(disturbedSurvival(robustType,totalTermType2Count,affectedTermType2Count) );
    }
    linRobust2 = mean(fitnessUponDeath);
    return linRobust2;
}

double Lineages::robustness1Gamma(int robustType,double shapeK, double meanU)
{
	if(linRobust1 > 0.0) return linRobust1;
    //if(internalTypeCnt <= 1) return 0.0;
    if(allLineages.size() < 1) coalescent();

    int sampleCnt = allLineages.size()*10;

    boost::mt19937 rng;
    boost::gamma_distribution<> gd(shapeK);
	boost::variate_generator<boost::mt19937&,boost::gamma_distribution<> > var_gamma(rng, gd);
	// meanU / shapeK * var_gamma()

    vector<int> tmpVec;
    map<int,int> affectedTermType2Count; // number of each types at terminal of a sub-lineage
    vector<double> fitnessUponDeath;
    map<int,int>::iterator it2;
    map<int,int> totalTermType2Count; // total number of each types at terminal
    vector<int> affectedDescendant;
    vector<int>::iterator it3;
    double rndDbl;
    for(unsigned int i=0; i<termTypes.size(); i++)
    {
        if(termTypes[i] != sparableCellType)
        {
            it2 = totalTermType2Count.find(termTypes[i]);
            if(it2 == totalTermType2Count.end())
            {
                totalTermType2Count[termTypes[i]] = 0;
            }
            totalTermType2Count[termTypes[i]] ++;
        }
    }
    if(totalTermType2Count.size() < 1) // none of the termini is important
    {
        return 0.0;
    }
    for(int i=0; i<sampleCnt; i++)
    {
        affectedDescendant.clear();
        affectedTermType2Count.clear();

        for(unsigned int j=0;j<(allIndex2DescendantIndex.size()-1);j++) // don't delete root
        {
			rndDbl = random(0,1);
			if(rndDbl < meanU / shapeK * var_gamma()) // delete this cell
			{
				tmpVec = allIndex2DescendantIndex[j];
				affectedDescendant.insert(affectedDescendant.end(),tmpVec.begin(),tmpVec.end());
			}
        }
		sort(affectedDescendant.begin(),affectedDescendant.end());
		it3 = unique(affectedDescendant.begin(),affectedDescendant.end());
		affectedDescendant.resize(it3 - affectedDescendant.begin());
		for(unsigned int j=0; j<affectedDescendant.size(); j++)
        {
            if(termTypes[affectedDescendant[j]] == sparableCellType) continue;
            it2 = affectedTermType2Count.find(termTypes[affectedDescendant[j]]);
            if(it2 == affectedTermType2Count.end())
            {
                affectedTermType2Count[termTypes[affectedDescendant[j]]] = 0;
            }
            affectedTermType2Count[termTypes[affectedDescendant[j]]] ++;
        }
        fitnessUponDeath.push_back(disturbedSurvival(robustType,totalTermType2Count,affectedTermType2Count) );
    }
	linRobust1 = mean(fitnessUponDeath);
	return linRobust1;
}

double Lineages::disturbedSurvival(int robustType, map<int,int> termType2Count, map<int,int> termType2DeathCount)
{
    double tmpDbl = 1.0;
    int tmpInt1,tmpInt2;
    map<int,int>::iterator it,it2;
    if(robustType == 1)
    {
        for(it=termType2Count.begin(); it!=termType2Count.end(); it++)
        {
            if(it->first == sparableCellType) continue;
            it2 = termType2DeathCount.find(it->first);
            if(it2 != termType2DeathCount.end())
            {
                tmpDbl = tmpDbl * (1.0 - ((double)it2->second / (double)it->second));
            }
        }
        return tmpDbl;
    }
    else if(robustType == 2)
    {
        tmpInt1 = 0;
        tmpInt2 = 0;
        for(it=termType2Count.begin(); it!=termType2Count.end(); it++)
        {
            if(it->first == sparableCellType) continue;
            tmpInt1 += it->second;
        }
        for(it=termType2DeathCount.begin(); it!=termType2DeathCount.end(); it++)
        {
            if(it->first == sparableCellType) continue;
            tmpInt2 += it->second;
        }
        return (1 - ((double)tmpInt2 / (double)tmpInt1));
    }
    else if(robustType == 3)
    {
        if(types2Weight.size() == 0)
        {
            cerr << "Define types2Weight before caculation of type 3 disturbedSurvival" << endl;
            exit(1);
        }
        for(it=termType2Count.begin(); it!=termType2Count.end(); it++)
        {
            if(it->first == sparableCellType) continue;
            it2 = termType2DeathCount.find(it->first);
            if(it2 != termType2DeathCount.end())
            {
                tmpDbl = tmpDbl * pow((1.0 - ((double)it2->second / (double)it->second)) , types2Weight[it->first]);
            }
        }
        return tmpDbl;
    }
    return 0.0;
}

double Lineages::robustness1(int robustType)
{
    if(termTypes.size() <= 1) return 0;
    //if(linRobust1 > 0.0) return linRobust1;
    if(allLineages.size() < 1) coalescent();

    //  map<int,vector<int> >::iterator it;
    map<int,int>::iterator it2;
    vector<int> tmpVec;
    map<int,int> affectedTermType2Count; // number of each types at terminal of a sub-lineage
    vector<double> fitnessUponDeath;

    for(unsigned int i=0; i<termTypes.size(); i++)
    {
        if(termTypes[i] != sparableCellType)
        {
            it2 = totalTermType2Count.find(termTypes[i]);
            if(it2 == totalTermType2Count.end())
            {
                totalTermType2Count[termTypes[i]] = 0;
            }
            totalTermType2Count[termTypes[i]] ++;
        }
    }
    if(totalTermType2Count.size() < 1) // none of the termini is important
    {
        return 0.0;
    }

    for(unsigned int i=0;i<(allIndex2DescendantIndex.size()-1);i++) // don't delete root
    {
        tmpVec = allIndex2DescendantIndex[i];
        affectedTermType2Count.clear();
        for(unsigned int i = 0; i<tmpVec.size(); i++)
        {
            if(termTypes[tmpVec[i]] == sparableCellType) continue;
            it2 = affectedTermType2Count.find(termTypes[tmpVec[i]]);
            if(it2 == affectedTermType2Count.end())
            {
                affectedTermType2Count[termTypes[tmpVec[i]]] = 0;
            }
            affectedTermType2Count[termTypes[tmpVec[i]]] ++;
        }
        fitnessUponDeath.push_back(disturbedSurvival(robustType,totalTermType2Count,affectedTermType2Count));
    }

    if(fitnessUponDeath.size() == 0) linRobust1 = 0.0;
    linRobust1 = mean(fitnessUponDeath);
    linRobustSE = stdev(fitnessUponDeath) / pow(fitnessUponDeath.size(),0.5);
    return linRobust1;
}

double Lineages::robustness1_Detail(int robustType, vector<double> &allRobustness)
{
    if(termTypes.size() <= 1) return 0;
    if(linRobust1 > 0.0) return linRobust1;
    if(allLineages.size() < 1) coalescent();

    //  map<int,vector<int> >::iterator it;
    map<int,int>::iterator it2;
    vector<int> tmpVec;
    map<int,int> affectedTermType2Count; // number of each types at terminal of a sub-lineage
    vector<double> fitnessUponDeath;

    for(unsigned int i=0; i<termTypes.size(); i++)
    {
        if(termTypes[i] != sparableCellType)
        {
            it2 = totalTermType2Count.find(termTypes[i]);
            if(it2 == totalTermType2Count.end())
            {
                totalTermType2Count[termTypes[i]] = 0;
            }
            totalTermType2Count[termTypes[i]] ++;
        }
    }
    if(totalTermType2Count.size() < 1) // none of the termini is important
    {
        return 0.0;
    }

    for(unsigned int i=0;i<(allIndex2DescendantIndex.size()-1);i++) // don't delete root
    {
        tmpVec = allIndex2DescendantIndex[i];
        affectedTermType2Count.clear();
        for(unsigned int i = 0; i<tmpVec.size(); i++)
        {
            if(termTypes[tmpVec[i]] == sparableCellType) continue;
            it2 = affectedTermType2Count.find(termTypes[tmpVec[i]]);
            if(it2 == affectedTermType2Count.end())
            {
                affectedTermType2Count[termTypes[tmpVec[i]]] = 0;
            }
            affectedTermType2Count[termTypes[tmpVec[i]]] ++;
        }
        fitnessUponDeath.push_back(disturbedSurvival(robustType,totalTermType2Count,affectedTermType2Count));
    }

    if(fitnessUponDeath.size() == 0) linRobust1 = 0.0;
    linRobust1 = mean(fitnessUponDeath);
    linRobustSE = stdev(fitnessUponDeath) / pow(fitnessUponDeath.size(),0.5);
    allRobustness = fitnessUponDeath;
    return linRobust1;
}

double Lineages::robustness2(int robustType)
{
    if(linRobust2 > 0.0) return linRobust2;
    if(internalTypeCnt <= 1) return 0.0;
    if(allLineages.size() < 1) coalescent();

    vector<int> allRulesOri;
    vector<int> allRules;
    for(int i=0; i<(internalTypeCnt-1); i++) allRulesOri.push_back(i);
    vector<int> pickedCells,tmpVec;

    map<int,int>::iterator it2;
    vector<int> affectedDescendant;
    map<int,int> affectedTermType2Count; // number of each types at terminal of a sub-lineage
    vector<double> fitnessUponDeath;

    for(unsigned int i=0; i<termTypes.size(); i++)
    {
        if(termTypes[i] != sparableCellType)
        {
            it2 = totalTermType2Count.find(termTypes[i]);
            if(it2 == totalTermType2Count.end())
            {
                totalTermType2Count[termTypes[i]] = 0;
            }
            totalTermType2Count[termTypes[i]] ++;
        }
    }
    if(totalTermType2Count.size() < 1) // none of the termini is important
    {
        return 0.0;
    }

    for(unsigned int i=0; i<allRulesOri.size(); i++)
    {
        affectedDescendant.clear();
        affectedTermType2Count.clear();
        pickedCells = internalTypes2allIndex[allRulesOri[i]];
        for(unsigned int j = 0; j<pickedCells.size(); j++)
        {
            tmpVec = allIndex2DescendantIndex[pickedCells[j]];
            affectedDescendant.insert(affectedDescendant.end(),tmpVec.begin(),tmpVec.end());
        }

        for(unsigned int j=0; j<affectedDescendant.size(); j++)
        {
            if(termTypes[affectedDescendant[j]] == sparableCellType) continue;
            it2 = affectedTermType2Count.find(termTypes[affectedDescendant[j]]);
            if(it2 == affectedTermType2Count.end())
            {
                affectedTermType2Count[termTypes[affectedDescendant[j]]] = 0;
            }
            affectedTermType2Count[termTypes[affectedDescendant[j]]] ++;
        }
        fitnessUponDeath.push_back(disturbedSurvival(robustType,totalTermType2Count,affectedTermType2Count) );
    }
    linRobust2 = mean(fitnessUponDeath);
    linRobustSE = stdev(fitnessUponDeath) / pow(fitnessUponDeath.size(),0.5);
    return linRobust2;
}

double Lineages::robustness2Sample(int robustType)
{
    if(linRobust2 > 0.0) return linRobust2;
    if(internalTypeCnt <= 1) return 0.0;
    if(allLineages.size() < 1) coalescent();
    if(ruleDeletionRate < 0)
    {
        cerr << "Define ruleDeletionRate first!" << endl;
        exit(1);
    }

    int sampleCnt = allLineages.size()*10;
    vector<int> allRulesOri;
    vector<int> allRules;
    for(int i=0; i<(internalTypeCnt-1); i++) allRulesOri.push_back(i);
    //
    // use (internalTypeCnt-1) so that you won't remove the rule governing the first division from the fertilized egg
    //

    vector<int> pickedRules;
    vector<int> pickedCells,tmpVec;
    double rndDbl;

    map<int,vector<int> >::iterator it;
    map<int,int>::iterator it2;
    map<int,int> totalTermType2Count; // total number of each types at terminal
    vector<int> affectedDescendant;
    vector<int>::iterator it3;
    map<int,int> affectedTermType2Count; // number of each types at terminal of a sub-lineage
    vector<double> fitnessUponDeath;
    int delCnt = 0;
    ptrdiff_t pd;

    for(unsigned int i=0; i<termTypes.size(); i++)
    {
        if(termTypes[i] != sparableCellType)
        {
            it2 = totalTermType2Count.find(termTypes[i]);
            if(it2 == totalTermType2Count.end())
            {
                totalTermType2Count[termTypes[i]] = 0;
            }
            totalTermType2Count[termTypes[i]] ++;
        }
    }
    if(totalTermType2Count.size() < 1) // none of the termini is important
    {
        return 0.0;
    }

    for(int i=0; i<sampleCnt; i++)
    {
        pickedRules.clear();
        affectedDescendant.clear();
        affectedTermType2Count.clear();
        rndDbl = random(0,1);
        for(unsigned int j=0; j<allRulesOri.size(); j++)
        {
            rndDbl -= binom(ruleDeletionRate,allRulesOri.size(),j);
            if(rndDbl < 0)
            {
                delCnt = j;
                break;
            }
        }

        allRules = allRulesOri;
        while(pickedRules.size() < (unsigned)delCnt)
        {
            pd = (ptrdiff_t)random(0,allRules.size());
            it3 = allRules.begin() + pd;
            pickedRules.push_back(*it3);
            allRules.erase(it3);
        }


        for(unsigned int j = 0; j<pickedRules.size(); j++)
        {
            pickedCells = internalTypes2allIndex[pickedRules[j]];
            for(unsigned int k = 0; k<pickedCells.size(); k++)
            {
                tmpVec = allIndex2DescendantIndex[pickedCells[k]];
                affectedDescendant.insert(affectedDescendant.end(),tmpVec.begin(),tmpVec.end());
            }
        }
        sort(affectedDescendant.begin(),affectedDescendant.end());
        it3 = unique(affectedDescendant.begin(),affectedDescendant.end());
        affectedDescendant.resize(it3 - affectedDescendant.begin());
        for(unsigned int j=0; j<affectedDescendant.size(); j++)
        {
            if(termTypes[affectedDescendant[j]] == sparableCellType) continue;
            it2 = affectedTermType2Count.find(termTypes[affectedDescendant[j]]);
            if(it2 == affectedTermType2Count.end())
            {
                affectedTermType2Count[termTypes[affectedDescendant[j]]] = 0;
            }
            affectedTermType2Count[termTypes[affectedDescendant[j]]] ++;
        }
        fitnessUponDeath.push_back(disturbedSurvival(robustType,totalTermType2Count,affectedTermType2Count) );
    }
    linRobust2 = mean(fitnessUponDeath);
    return linRobust2;
}

double Lineages::robustness2Sample_Detail(int robustType, vector<double> & allRobustness)
{
    if(linRobust2 > 0.0) return linRobust2;
    if(internalTypeCnt <= 1) return 0.0;
    if(allLineages.size() < 1) coalescent();
    if(ruleDeletionRate < 0)
    {
        cerr << "Define ruleDeletionRate first!" << endl;
        exit(1);
    }

    int sampleCnt = allLineages.size()*10;
    vector<int> allRulesOri;
    vector<int> allRules;
    for(int i=0; i<(internalTypeCnt-1); i++) allRulesOri.push_back(i);
    //
    // use (internalTypeCnt-1) so that you won't remove the rule governing the first division from the fertilized egg
    //

    vector<int> pickedRules;
    vector<int> pickedCells,tmpVec;
    double rndDbl;

    map<int,vector<int> >::iterator it;
    map<int,int>::iterator it2;
    map<int,int> totalTermType2Count; // total number of each types at terminal
    vector<int> affectedDescendant;
    vector<int>::iterator it3;
    map<int,int> affectedTermType2Count; // number of each types at terminal of a sub-lineage
    vector<double> fitnessUponDeath;
    int delCnt = 0;
    ptrdiff_t pd;

    for(unsigned int i=0; i<termTypes.size(); i++)
    {
        if(termTypes[i] != sparableCellType)
        {
            it2 = totalTermType2Count.find(termTypes[i]);
            if(it2 == totalTermType2Count.end())
            {
                totalTermType2Count[termTypes[i]] = 0;
            }
            totalTermType2Count[termTypes[i]] ++;
        }
    }
    if(totalTermType2Count.size() < 1) // none of the termini is important
    {
        return 0.0;
    }

    for(int i=0; i<sampleCnt; i++)
    {
        pickedRules.clear();
        affectedDescendant.clear();
        affectedTermType2Count.clear();
        rndDbl = random(0,1);
        for(unsigned int j=0; j<allRulesOri.size(); j++)
        {
            rndDbl -= binom(ruleDeletionRate,allRulesOri.size(),j);
            if(rndDbl < 0)
            {
                delCnt = j;
                break;
            }
        }

        allRules = allRulesOri;
        while(pickedRules.size() < (unsigned)delCnt)
        {
            pd = (ptrdiff_t)random(0,allRules.size());
            it3 = allRules.begin() + pd;
            pickedRules.push_back(*it3);
            allRules.erase(it3);
        }


        for(unsigned int j = 0; j<pickedRules.size(); j++)
        {
            pickedCells = internalTypes2allIndex[pickedRules[j]];
            for(unsigned int k = 0; k<pickedCells.size(); k++)
            {
                tmpVec = allIndex2DescendantIndex[pickedCells[k]];
                affectedDescendant.insert(affectedDescendant.end(),tmpVec.begin(),tmpVec.end());
            }
        }
        sort(affectedDescendant.begin(),affectedDescendant.end());
        it3 = unique(affectedDescendant.begin(),affectedDescendant.end());
        affectedDescendant.resize(it3 - affectedDescendant.begin());
        for(unsigned int j=0; j<affectedDescendant.size(); j++)
        {
            if(termTypes[affectedDescendant[j]] == sparableCellType) continue;
            it2 = affectedTermType2Count.find(termTypes[affectedDescendant[j]]);
            if(it2 == affectedTermType2Count.end())
            {
                affectedTermType2Count[termTypes[affectedDescendant[j]]] = 0;
            }
            affectedTermType2Count[termTypes[affectedDescendant[j]]] ++;
        }
        fitnessUponDeath.push_back(disturbedSurvival(robustType,totalTermType2Count,affectedTermType2Count) );
    }
    linRobust2 = mean(fitnessUponDeath);
    allRobustness = fitnessUponDeath;
    return linRobust2;
}

void Lineages::coalescent()
{
    internalTypeCnt = 0;
    vector<string> type2desc; // map internal cell type to description string
    map<string,int> lineage2type; // map lineage string to internal/terminal cell type
    map<string,int>::iterator it,it2;
    int currentIndex = 0;

    char buffer [50]; // for conversion of number->string
    bool newType = false;
    string tmpDesc,tmpSibling1,tmpSibling2,tmpParent,parentDesc1,parentDesc2;
    //int testInt;
    //string testStr;
    vector<int> testVec;

    //populate the mapping
    for(unsigned int i=0; i<termTypes.size(); i++)
    {
        //testStr = termLineages[i];
        lineage2type[termLineages[i]] = termTypes[i];
        lineages2allIndex[termLineages[i]] = i;
        allLineages.push_back(termLineages[i]);
        allIndex2DescendantIndex[i].push_back(i);
    }

    while(lineage2type.size()>1)
    {
        //find siblings for a random node
        for ( it=lineage2type.begin() ; it != lineage2type.end(); it++ )
        {
            tmpSibling1 = it->first;
            tmpParent = tmpSibling1.substr(0,int(tmpSibling1.length()-1));
            tmpSibling2 = tmpSibling1.substr(int(tmpParent.length()),1);
            if(tmpSibling2 == "1")
            {
                tmpSibling2 = "0";
                tmpSibling2.insert(0,tmpParent);
            }
            else
            {
                tmpSibling2 = "1";
                tmpSibling2.insert(0,tmpParent);
            }
            it2 = lineage2type.find(tmpSibling2);
            if(it2 == lineage2type.end()) continue; // could not find the other sibling
            currentIndex = allLineages.size();
            allLineages.push_back(tmpParent); // found the other sibling
            lineages2allIndex[tmpParent] = currentIndex;

            //testVec = allIndex2DescendantIndex[currentIndex];
            //testInt = lineages2allIndex[tmpSibling1];
            allIndex2DescendantIndex[currentIndex].insert
            (
                allIndex2DescendantIndex[currentIndex].end(),
                allIndex2DescendantIndex[lineages2allIndex[tmpSibling1]].begin(),
                allIndex2DescendantIndex[lineages2allIndex[tmpSibling1]].end()
            );
            allIndex2DescendantIndex[currentIndex].insert
            (
                allIndex2DescendantIndex[currentIndex].end(),
                allIndex2DescendantIndex[lineages2allIndex[tmpSibling2]].begin(),
                allIndex2DescendantIndex[lineages2allIndex[tmpSibling2]].end()
            );
            //testVec = allIndex2DescendantIndex[currentIndex];
            break;
        }

        //determine the parent type, record if necessary
        newType = true;

        parentDesc1 = "";
        sprintf(buffer,"%d",lineage2type[tmpSibling1]);
        parentDesc1.insert(0,buffer);
        parentDesc1.insert(0,"_");
        sprintf(buffer,"%d",lineage2type[tmpSibling2]);
        parentDesc1.insert(0,buffer);

        parentDesc2 = "";
        sprintf(buffer,"%d",lineage2type[tmpSibling2]);
        parentDesc2.insert(0,buffer);
        parentDesc2.insert(0,"_");
        sprintf(buffer,"%d",lineage2type[tmpSibling1]);
        parentDesc2.insert(0,buffer);

        for(unsigned int i=0; i<type2desc.size(); i++)
        {
            if((type2desc[i] == parentDesc1) or (type2desc[i] == parentDesc2)) //if not new type
            {
                newType = false;
                internalTypes.push_back(i);
                internalDesc.push_back(parentDesc1);

                if(lineage2type[tmpSibling1] == lineage2type[tmpSibling2]) {
                    prolifDepth.push_back(tmpParent.length());
                } else {
                    diffDepth.push_back(tmpParent.length());
                }

                //combine them into one single node
                lineage2type.erase(tmpSibling1);
                lineage2type.erase(tmpSibling2);
                lineage2type[tmpParent] = i;
                internalTypes2allIndex[i].push_back(currentIndex);

                break;
            }
        }
        if(newType)
        {
            //testInt = type2desc.size();
            //testInt = internalTypes.size();
            type2desc.push_back(parentDesc1);
            //testStr = type2desc[0];
            //testInt = type2desc.size();
            internalDesc.push_back(parentDesc1);
            internalTypes.push_back(internalTypeCnt);

            if(lineage2type[tmpSibling1] == lineage2type[tmpSibling2]) {
                prolifDepth.push_back(tmpParent.length());
            } else {
                diffDepth.push_back(tmpParent.length());
            }

            //combine them into one single node
            lineage2type.erase(tmpSibling1);
            lineage2type.erase(tmpSibling2);
            lineage2type[tmpParent] = internalTypeCnt;
            internalTypes2allIndex[internalTypeCnt].push_back(currentIndex);

            internalTypeCnt ++;


        }
        currentIndex ++;
    }//repeat until only one node left

    //testInt = termTypes.size();
    //testInt = type2desc.size();
    //testInt = internalTypes.size();
    testVec = internalTypes2allIndex[0];
    testVec = internalTypes2allIndex[1];
    testVec = internalTypes2allIndex[2];
    linCompl = (((double)type2desc.size()) / ((double)internalTypes.size()) );
}

void Lineages::getDepthStat()
{
    if(allLineages.size() < 1) coalescent();

    double tmpDepth;
    vector<int> allTermDepth;
    for(unsigned int i=0; i<termLineages.size(); i++)
    {
        tmpDepth = termLineages[i].length();
        allTermDepth.push_back((int)tmpDepth);
        if(tmpDepth  > depthMax) depthMax = tmpDepth;
    }
    depthMean =  mean(allTermDepth);
    depthSD = stdev(allTermDepth);

    vector<int> depth2InternCnt;
    vector<int> depth2TermCnt;
    for(int i=0; i<=depthMax; i++)
    {
        depth2InternCnt.push_back(0);
        depth2TermCnt.push_back(0);
    }
    for(unsigned int i=0; i<allLineages.size(); i++)
    {
        if(allIndex2DescendantIndex[lineages2allIndex[allLineages[i]]].size() == 1) continue;
        depth2InternCnt[allLineages[i].length()] ++;
    }
    for(unsigned int i=0; i<termLineages.size(); i++)
    {
        depth2TermCnt[termLineages[i].length()] ++;
    }
    depth2DivRate.clear();
    for(int i=0; i<=depthMax; i++)
    {
        depth2DivRate.push_back((double)depth2InternCnt[i]/(double)(depth2TermCnt[i] + depth2InternCnt[i]));
    }
    return;
}

double Lineages::factln(int n)
{
    double x = (double)n + 1.0;
    double tmp = x + 5.5;
    tmp -= (x + 0.5) * log(tmp);
    double ser = 1.000000000190015
                 + 76.18009172947146    / (x + 1)
                 - 86.50532032941677    / (x + 2)
                 + 24.01409824083091    / (x + 3)
                 -  1.231739572450155   / (x + 4)
                 +  0.12086509738661e-2 / (x + 5)
                 -  0.5395239384953e-5  / (x + 6);
    return log(2.5066282746310005 * ser/(x) ) - tmp;
}

double Lineages::binom(double p, int n, int k)
{
    return exp(factln(n) - factln(k) - factln(n - k) + k * log(p) + (n-k) * log(1-p));
}
